# This script generates AR data corresponding to a model.
# It requires the package "mAr" located on CRAN - run install.packages("mAr")
# See http://www.oga-lab.net/RGM2/functions.php?show=all&query=package:mAr
#
# Form the autoregressive coefficients and check for stationarity. 
# There are two conditions which we are trying to meet:
#
# 1. sparsity in the spectral density matrix (ie. autocorrelations between the variables in the
# time series) which reduces to sparsity patterns in products of the coefficients - see Dalhaus
# 2000, pp. 164. 
#
# 2. Weak stationarity. The stationarity condition is equivalent to eigenvalues of the companion matrix
# being outside the unit circle. See faculty.washington.edu/ezivot/econ584/notes/varModels.pdf
# for the details (pg. 385). 
#
# If there is a better way to ensure stationarity (perhaps by generating in the spectral domain?)
# please let me know. Perhaps there is a stationarity condition which relies solely on the spectral 
# density matrix, so that we can start by assuming a spectral density matrix which, for example, is 
# banded, and then use the Cholesky factors to recover the AR coefficients? 
###################################################################################################

source('check_sparsity.R')
library(mAr)

# p is the data dimension, N is the length of the time series, C is the error covariances,
# w is the intercept term (0 for us), and A is a p x pk dimensional matrix where k is the lag
# with A = [A_1,...,A_p]. The sequence generated will be x[t] = \sum_{i} A_i X[t-i]. 
# Note that is the opposite sign convention as that used in Songsiri et al. 
p   <- 3
N   <- 10000
C   <- diag(rep(1, p))
w   <- rep(0, p)
A_1 <- diag(rep(0.5, p))
A_2 <- diag(rep(0.3, p))
A_3 <- diag(rep(0.1, p))
A   <- cbind(A_1, A_2, A_3)

# Check the sparsity structure of the inverse spectral density matrix. 

temp <- check_sparsity(A, C, p = p, k = 3)

# Form companion  matrix F and check for stationarity. F will be pk x pk. 
# The stationarity condition is that the eigenvalues of F are within the complex
# unit ball. 
F <- matrix(0, 3 * p, 3 * p)
F[1:p, 1:(3 * p)] <- A
diag(F[-(1:p), ]) <- 1
f.eig <- eigen(F)$values
stationary.check <- ( max(abs(f.eig)) < 1 )


# Generate data using mAr.sim

Data <- mAr.sim(w, A, C, N)


###################################################################################################
# This is a more complicated example of how to choose a sparsity structure, then enforce stationarity. 
# We follow the procedure in Songsiri et al., Experiment 2, pg. 2689.


